In [1]:
%matplotlib notebook
from astropy.cosmology import LambdaCDM
import numpy as np
import matplotlib.pyplot as plt
from astropy import constants as const
import astropy.units as u
from scipy.integrate import quad
In [2]:
cosmo = LambdaCDM(H0=70, Om0=0.3, Ode0=0.7, Tcmb0=2.725)
In [3]:
z = np.arange(0, 1.5, 0.01)
z=0.4
In [4]:
phi_star = 3.6 * cosmo.efunc(z)**2
alpha = -1.05 * (1 + z)**(-2/3)
fr = 0.8*(1 + z)**(-1/2)
In [5]:
alpha
Out[5]:
Need to compute the cluster volume...
$M_{vir} = 4/3 \pi r^3_{vir} \rho_c(r<r_{vir}) = 4/3 \pi r^3_{vir} \Delta_c \rho_c$
if we let $\Delta_c = 200$ then
$M_{200} = 4/3 \pi r^3_{200} 200 \rho_c$ with $\rho_c = \frac{3H(z)^2}{8\pi G}$
or just $M_{200} = V_{200}200\rho_c$
Don't forget that $H(z) = H_0E(z)$
In [6]:
def rho_crit(z, cosmo):
# convert G into better units:
G = const.G.to(u.km**2 * u.Mpc/(u.M_sun * u.s**2))
return 3 / (8 * np.pi * G) * cosmo.H0**2 * cosmo.efunc(z)**2 # Mpc^3
# So now we are going to calculate the volumes as a function of z
M200 = 1e15 * u.solMass
V200 = M200/ (200 * rho_crit(z, cosmo))
In [7]:
V200
Out[7]:
The Schechter Function:
For Luminosity:
$\Phi(L) = \phi^\star \frac{L}{L_\star}^\alpha e^{-\frac{L}{L_\star}}$
For Magnitudes:
$\Phi(M) = \phi^\star\frac{2}{5}log(10) (10^{\frac{2}{5}(M_\star - M)})^{\alpha+1} e^{-10^{\frac{2}{5}(M_\star - M)}}$
In [8]:
def schechterL(luminosity, phiStar, alpha, LStar):
"""Schechter luminosity function."""
LOverLStar = (luminosity/LStar)
return (phiStar/LStar) * LOverLStar**alpha * np.exp(- LOverLStar)
def schechterM(magnitude, phiStar, alpha, MStar):
"""Schechter luminosity function by magnitudes."""
MStarMinM = 0.4 * (MStar - magnitude)
try:
return (0.4 * np.log(10) * phiStar * 10.0**(MStarMinM * (alpha + 1.)) * np.exp(-10.**MStarMinM))
except OverflowError:
return 0.0
In [9]:
M200 = 10 # 10^14 Msun
Mpiv = 6 # 10^14
zpiv = 0.6
alpha = -0.96 * (M200 / Mpiv)**0.01 * ((1 + z)/ (1 + zpiv))**-0.94
Phi = 1.68 * (M200 / Mpiv)**0.09 * ((1 + z)/ (1 + zpiv))**0.09 * cosmo.efunc(z)**2
fr = 0.62 * (M200 / Mpiv)**0.08 * ((1 + z)/ (1 + zpiv))** -0.80
In [10]:
M_star = -20.44 + 5 * np.log10(0.7) # abs mag.
M_star = -21.946521706412092 # abs mag @ z= 0.4
M_star_sub = M_star - 2.5 * np.log10(0.4)
M_star_big = M_star - 2.5 * np.log10(100)
#M_star = 20.071802735723015 # app mag @ z= 0.4
y, err = quad(schechterM, -np.inf, -21.61828907152557, args=(phi_star, alpha, M_star))
#y, err = quad(schechterM, -np.inf, M_star_sub, args=(phi_star, alpha, M_star))
In [11]:
(y * V200.value + 1) * fr
Out[11]:
In [12]:
Phi
Out[12]:
In [13]:
phi_star
Out[13]:
In [14]:
1/5
Out[14]:
Here is the stuff from the c++ code.
In [22]:
def dAintegrand2(y, alpha, mstar):
#val = y**alpha * np.exp(-y)
# using mags instead of luminosity
val = 10.0**(-2 / 5 * (y - mstar) * (alpha + 1.)) * np.exp(-10.** (-2 / 5) * (y - mstar))
return val
def calc_hon(z, mlimit, mstar):
M200 = 10 # 10^14 Msun
Mpiv = 6 # 10^14
zpiv = 0.6
alpha = -0.96 * (M200 / Mpiv)**0.01 * ((1 + z)/ (1 + zpiv))**-0.94
Phi = 1.68 * (M200 / Mpiv)**0.09 * ((1 + z)/ (1 + zpiv))**0.09 * cosmo.efunc(z)**2
fr = 0.62 * (M200 / Mpiv)**0.08 * ((1 + z)/ (1 + zpiv))** -0.80
R200 = M200 * 4301.8644383/ (cosmo.H0.value**2 * cosmo.efunc(z)**2)
R200 = R200**(1 / 3)
# Calculate factors
Vol = 4 / 3 * np.pi * R200**3 # This gives basically the same number as my calculation above.
# Using Mags Instead of Luminosity
NN = -0.4 * np.log(10) * Vol * Phi * quad(dAintegrand2, -np.inf, mlimit, args=(alpha, mstar))[0]
print(NN)
return (NN + 1) * fr
In [23]:
calc_hon(z, -21, M_star)
Out[23]:
In [17]:
type(np.log(10))
Out[17]:
In [18]:
cosmo.efunc(z)**2
Out[18]:
In [19]:
from scipy import special
In [20]:
phi_star * special.gammainc(alpha + 1, 3) * V200 * fr
Out[20]:
In [21]:
alpha
Out[21]:
In [ ]: